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We investigate a scenario in which neutrinos are coupled to a pseudoscalar degree of freedom 
and where decays v\ — \ v% + ip and inverse decays are the responsible mechanism for obtain- 
ing equilibrium. In this context we discuss the implication of the invisible neutrino decay on the 
neutrino-pseudoscalar coupling constant and the neutrino lifetime. Assuming the realistic scenario 
of a thermal background of neutrinos and pseudoscalar we update the bound on the (off-diagonal) 
neutrino-pseudoscalar coupling constant to g < 2.6 x 10~ 13 and the bound on the neutrino life- 
time to r < 1 x 10 13 s. Furthermore we confirm analytically that kinetic equilibrium is delayed by 
(2jT)i two Lorentz 7-factors - one for time dilation of the (decaying) neutrino lifetime and one from the 

^"t ' opening angle. We have also confirmed this behavior numerically. 
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O ■ I. INTRODUCTION 

u 

The possibility for neutrino interactions beyond the standard model has been studied in many contexts over the 
years. One particularly simple possibility is that neutrinos couple to a new pseudoscalar degree of freedom, as is, for 
1 example, the case in Majoron models - see the following references for previous discussions about the dynamics of the 
strong neutrino-pseudoscalar coupling and its astrophysical implications [1-10]. 

Astrophysics provides fairly stringent constraints on such couplings. For example SN1987A provides a bound on 
the dimensionless coupling constant of order 10 -7 < g < 10~ 5 [11-14] by requiring that the neutrino signal should not 
— — be significantly shortened. 

In the same way there are two cosmological bounds on g. First, the value of g should not be large enough that 
pseudoscalars are fully thermalized before big bang nucleosynthesis. This leads to g < 10 -5 . Second, a significant 
value of g will make neutrinos self-interacting in the late universe and prevent neutrino free-streaming. This possibility 
I has been discussed a number of times in the literature (see [15-19]). 

The effect on cosmological observables such as the cosmic microwave background (CMB) spectrum were studied 
in [15-17, 19-23]); particularly it was found that although models with no neutrino free-streaming can mimic the 
__, . matter power spectrum of ACDM models, they produce a distinct signature in the CMB spectrum which is much 
t— ( ' harder to reproduce. The feature arises because neutrinos act as a source term for photon perturbations. If there is 
0^ , no free-streaming, the source term is stronger and consequently the CMB anisotropy is increased for all scales inside 
' the particle horizon at recombination. On the other hand, there is no effect on larger scales. 

This distinct signature has been used to constrain models without neutrino free-streaming and in [20, 23] it was 
used to constrain the corresponding neutrino-pseudoscalar coupling parameters. 

However, for decays and inverse decays the interaction was treated in a somewhat simplified manner in the sense 
that the momentum equilibration rate was assumed to be roughly T* ~ l/(7 2 r), where t is the rest- frame lifetime 
and 7 ~ E v /m v is the Lorentz boost factor. 

In this paper we wish to check this assumption in an explicit way with a realistic setup. The possible departure from 
the simple relation T* ~ 1/(7 2 t) is something which is highly relevant for parameter estimations - such as placing 
bounds on the neutrino-pseudoscalar coupling and hence also for constraining the neutrino lifetime. Furthermore it 
is something which needs to be taken into account in numerical studies in which we allow for nonstandard neutrino 
interactions. One particular area where detailed knowledge of the interaction would be useful is the search for the 
cosmic neutrino background [24-27]. 

One comment is in order here: The motivation for looking at decays and inverse decays rather than various scattering 
processes involving neutrinos and pseudoscalars (yv — > (pip, iptp — > vv, vip — > v<p) is that the probabilities of these 
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scattering processes are proportional to g 4 , where g is the neutrino-pseudoscalar coupling constant. The probability 
of the decay V\ — > v 2 + <p, on the other hand, is proportional to g 2 . Consequently, at small values of g the decay 
actually dominates over the scattering processes and allows us to put severe constraints on g. 

This paper is organized as follows: In Sec. II we look at the setup with a gas consisting of two neutrino species and 
a pseudoscalar - a gas that only has decays and inverse decays to obtain equilibrium. We argue for the — 7 in the 
decay rate. In Sec. Ill we look at an initial situation of a standing wave of the heavy neutrinos and no other particles. 
In Sec. IV we introduce thermal distributions of the light neutrino and of pseudoscalars into a thermal background 
while keeping the initial conditions for the heavy neutrino. Furthermore we discuss the implication of the decay on 
the neutrino-pseudoscalar coupling and on the neutrino lifetime. We present numerical results in Sec. V. Finally we 
have a conclusion and an appendix concerning calculations for the numerical implementation of the system. 



II. THERM ALIZATION OF A GAS WITH ONLY DECAYS AND INVERSE DECAYS 

Thermalization of a gas by decay and inverse decay is a nontrivial process because of phase space limitation. As 
long as one of the involved particles can interact with an external heat bath, it is in principle possible to thermalise 
the gas provided that the interaction rate is sufficiently fast. This is, for instance, the case with thermal leptogenesis 
in which the decay products are thermalized by SM gauge interactions. 

However, for the case studied here this is not true. The weak interactions are far too weak to maintain equilibrium 
at the eV temperatures considered here. In this case full thermal equilibrium can never be achieved. 

The standard case usually studied, for example, in the case of thermal leptogenesis is a spatially homogeneous gas 
in which interactions drive the distribution toward thermal equilibrium (see e.g. [28, 29]) 

However, from the point of view of structure formation and more specifically free-streaming the important point 
is the rate of directional momentum transfer between species. For example, Thompson scattering is inefficient for 
maintaining energy equilibration between electrons and photons, but very efficient for exchanging momentum between 
the two species. This can be seen from the simple relations \AE 7 /E 7 \ <~ E 1 /m e and \Ap/p\ <~ I in a single scattering 
event. Therefore Thompson scattering is very efficient for driving the acoustic baryon-photon oscillations prior to 
recombination. 

However, for a gas with only decays and inverse decays momentum transfer is even more inefficient than energy 
transfer. Roughly the energy transfer time scale is given by the decay rate T = 1/(71"), i.e. the usual Lorentz 
suppressed rest-frame decay rate. However, in the lab frame the decay products are emitted in a cone of opening 
angle I/7 relative to the direction of momentum of the parent particle. Therefore, in a single decay the momentum 
direction is changed by only \Ap/p\ ~ I/7. This finally means that the rate of momentum change in the gas is roughly 
1/(7 2 t); i.e. for relativistic decays it is highly suppressed, and even suppressed relative to the energy exchange rate. 

Let us begin with the Lagrangian for a generic pseudoscalar neutrino interaction 1 



C = -i^2g 3 k<fVjj5Vk- (1) 

j,k 

We will consider only two neutrinos, one we consider to be heavy (y\ [or just 1 for convenience] with mass mi), 
a massless neutrino [y 2 [or just 2 for convenience]) and a massless pseudoscalar (ip). Thus we drop the index of g 
(g = #1,2) and the Lagrangian becomes 

C = -igip{Vu 5 u 2 + (2) 

In the following section we derive the specific Boltzmann collision integrals relevant for decays and inverse decays 
in an inhomogeneous gas. 

The variation of any overall quantity Q can be calculated from the distribution functions: 



Voto = £y (2K?m (2ir)323, (2,)32^ (27r) * ^ ~ ^ ~ 

[f2.Ua - fi) - A(l - /2)(1 + U)]QiSi (3) 



1 The coupling structure could in principle be derivative instead of pseudoscalar. However, this point makes no difference to the discussion 
here since we study only decays and inverse decays. One could also choose a scalar coupling - it would only lead to a very small difference, 
which, in fact, is removed completely in the approximation where the lighter neutrino is massless. 
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where Si = 1 and S 2 = S v = — 1. 



III. TOTAL TRANSVERSE MOMENTUM - WITH NO BACKGROUND 

We want to calculate the initial transverse momentum when we start with a standing wave of l's. The distribution 
functions are 

h = y(« 3 (p1-p6) + * 3 (p1+po)) 

h = U = o. (4) 

where po is the momentum of the standing wave (and Eq will be the corresponding on-shell energy) . Since the two 
terms initially contribute equally, we are free to change to one beam instead 2 

/1 = ni<5 3 0! -po) 

h = U = 0- (5) 

This will give the same result. 

We calculate the matrix element. From tracing and averaging over incoming and summing over outgoing spins and 
assuming the masses of the neutrinos mi = m,m 2 = 0: 3 



So we find 



|M| 2 = 2 5 2 (pi-p 2 -mim 2 )=5 2 TO 2 . (6) 
(t) V- f d* Pl d 3 p 2 d 3 p. 



Ei u pi u p 2 a p v 4 4 

g 2 m 2 niS 3 ( Pl - po)QiSi 
-g 2 m 2 n 1 ^ f d 3 p 1 d 3 p 2 d 3 p ip 4 3 ^ 

What we want to find is a measure of the transverse momentum created in the very beginning. Obviously, there is 
no momentum if we just sum over the transverse momentum vectors, so we sum the magnitudes of created transverse 
momenta instead. This means putting —QiSi = \pl x p" \. Thus we have 

5 2 m 2 ni f d 3 p 1 d 3 p 2 d 3 p v> 4 3 ^ 

voto = "w J e iE2 e v 5 {P1 - P2 p - )s (P1 - Po) 

(\p~i x Pol + \P~2 x Po| + \p~tp x Pol) • (7) 

After integration we find 

d\P±\ 2 3 

dt g m n\ 



Volume 64(27r) 4 J B ' 
Inserting g 2 m = where t is the 1-lifetime, and = ^ our final result is 

d\P±\ „ j? , 



Volume 8(27r) 3 T7 2 tj 

as expected. 



<* 37^ (8) 



2 In the numerical calculation we do not change to one beam 

3 Had one chosen a scalar interaction, the mit7i2 term would change sign. However, since we set m,2 = anyway, it makes no difference 
under this approximation. 
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IV. A MORE GENERAL CASE 

In case of a thermal background of neutrinos and pseudoscalars we cannot use the simple approach specified by 
Eqs. 4 and 5. Hence we consider all the distribution functions 

- h) - /i(i - / 2 )(i - /„) = hU - A(i + - h). (9) 

This complicates things a bit. But if we again take only the initial time it is solvable. The term we are discussing is 
\f2f1p — + ftp — f2)\\Pi x Po|- We assume initial equilibrium densities of 2 and ip, which using Boltzmann statistics 

J^J Irp 

means that the distribution functions of 2 and <p are identical. It also means that we can use f 2 f v = c . Hence 

we can split the term in an f\ part and an extra part. The /1 part yields exactly the same result as in Sec. Ill, but 
the extra part is proportional to fif v and yields the following: 

f d 3 px d 3 p 2 d 3 p v 4 4 
- ' -(27r) 4 (5 4 (pi -p 2 -p v ) 



Volumeoxtra J (2tt) 3 2E 1 {2ir) 3 2E 2 {2ir) 3 2E v 
-E\/T 



9 2 m 2 c 1 (\p[ x p%\ + {pi x ph\ + \p v x pol) (10) 



which reduces to 



with relativistic limit 



9WA 



dt 



Volume extra 



d\P±\ 

dt 



(11) 



EqPq 



(E T + 2T 2 ) . (12) 



Thus the final result is 

9\P±\ 



dt _ 



Volume extra 

E oni , E oPo f^^-m/T f°° ^ -E./T 



Volume 8(2vr) 3 T7 2 
with relativistic limit 



E po ( -m/T f°° -E 1 /T \ , 



Volume 8(27r)3r 7 2 4^ 2 r 7 1 ° >' K ' 

Hence, we are led to conclude that in the realistic scenario of having a background thermal distribution of light 
neutrinos as well as of pseudoscalars, we see a correction in the form of the second term in Eq. 14. This extra 
contribution is something which should be taken into account when putting bounds on the neutrino-pseudoscalar 
interaction. 

We can make a rough estimate of the improvement on the bound of the decay coupling constant when taking Eq. 14 

into account. First, we notice that in the presence of the first term on the right-hand side of Eq. 14 we are investigating 

2 

the standard case which was previously studied in [20]. Here the naive decay rate Tdecay = is translated into a 

transport rate Transport = Tdccay (^) 3 ; where the factor (^) 3 is due to three Lorentz gamma factors: The first one 
comes from transforming from the rest frame of the parent neutrino to the frame of the thermal medium. This will 
give us a decay rate in the frame of the thermal medium. The other two come from the following reason: The decay 
is isotropic in the rest frame of the parent neutrino; however, the decay products will have directions within an angle 
corresponding to a factor 7. So, to randomize the direction of the original neutrino we must include another factor of 
7. In total when we transform from the medium frame decay rate to the relevant transport rate we get two factors of 
gamma. All in all we arrive at the desired expression 

/ TTl \ 

J- transport — rdccay^J • (15) 
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To ensure that the neutrinos are still free-streaming at the time of photon decoupling as required by observations of 
the CMB [15, 23], we can then compare the transport rate with the expansion rate of the universe -f/dec at photon 
decoupling. The requirement for free-streaming is T < Hdec- This leads to the bound [20] 

^o.eixio- 11 ^-^) 2 . (16) 

In the event of the decay taking place in a thermal distribution of light neutrinos and pseudoscalars we need to take 
the second term in Eq. 14 into account. Especially since at the time of photon decoupling, assuming a generic heavy 
neutrino mass of m v = 50meV and energy E = ?>T V ^ CC <~ 3 x 0.18 eV, we get 7 ~ ^ <~ 3.6. Combined with the fact 
that for a relativistic species {m v < 3T„ i dcc) and for our relativistic heavy neutrino n v ~ X^ dcc up to factors of order 
unity, this means that the transport rate we should be comparing is rather 



Transport ~ Tdccay ^(^g) ^^ 7F ("g) ^ ' 



(17) 



where the factor of 167r takes into account this missing factor in the denominator of the second term. If we translate 
into a bound on the decay coupling constant, this gives 



n n ,„ is /50meV\ 3/2 / T v \ 1/2 , /50meV\ / T v 

5 <2.6xl0~ 13 - " 1 + I.8xl0~ 3 ) ( TTTTTTT 

y \ m v J \0.18eV J \ \ m v ) V0.18eV 



-1/2 



For 7 <~ 3.6 with a neutrino mass m v = 50meV this translates simply into the bound 

5 <2.6xl0~ 13 , (19) 

i.e. an improvement of more than a factor of 10. Translating this into a limit on the neutrino lifetime in the restframe 
we get 

r<lxl0 13 s, (20) 

hence there is still the possibility for the neutrino to be short lived when we let the decay take place in a thermal 
background. 



V. NUMERICAL RESULTS 



In the numerical implementation we had to change the setup a little. Because of problems with the bins, it was 
impossible to get the heavy neutrino to have vanishing momentum in the transverse direction. This could not be 
remedied by increasing the number of bins. Therefore we made a thermal distribution of l's around an average 
momentum. Specifically we chose a standard scenario with m = 2T, (p x ) = ±3T and made a thermal distribution 
around this. We chose the artificially high value g = -^7= (making T = ^ 2L ) to let the code find equilibrium in a 
reasonable time. We checked that the code did reach equilibrium both under the initial condition of no light neutrinos 
from the beginning and under the assumption of an initial thermal distribution of light neutrinos. More on the 
numerics is provided in the Appendix VII. 

This means that the calculated formulas cannot be verified explicitly, since we have no well-defined gamma factor. 
However, the factor of ^ can almost be found. First, the ^ is, of course, trivial. If the coupling is weakened, the 

lifetime is correspondingly longer. We checked that our code yielded this result. The ^ for the time dilation of the 
neutrino lifetime is also quite trivial. We found this as well - by noticing the decreased numbers of 2's produced in 
the very first step when we used an alternative scenario: m = 2T, (p x ) = ±4T. Since we had thermal distributions 
rather than sharply defined momenta, we could not expect to find the exact relation between average 7 s to be the 
same as the relation between the created particles. 

However, the most interesting second i can be illustrated by numerical plots. Fig. 1 shows the first distributions 
after the very first step, whereas Fig. 2 shows the distributions after the first step in the alternative scenario with 
a larger gamma factor. Two things are very important to notice. First, the created 2's are indeed anisotropic. 
This is the effect of taking their isotropic distribution in the frame of the decaying particle 1 and making a Lorentz 
transformation to the (cosmic) laboratory frame. Second, the fact that the alternative scenario shows more anisotropy 
among the 2's should make it clear that it must be a gamma factor. One could alter Eqs. (13), (14) to match the 
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-50 50 -50 50 

px px 



FIG. 1: The situation just after the first step, FIG. 2: The situation just after the first step, 

t = 9.95 x 1CP 9 r, in the standard scenario. The lines are t = 9.95 x 1CP 9 r, in the alternative scenario. The lines are 
contour plots of the distribution function log(/ * T) - as contour plots of the distribution function log(/ * T) - as 

defined in the Appendix. p y means transverse momentum and defined in the Appendix. p y means transverse momentum and 
is measured in units of bin length dp = . is measured in units of bin length dp — j. 




initial conditions for the numerics. But since Eqs. (13), (14) are complicated enough already, and since the numerics 
is not done to create new results but only to confirm the pattern of Eqs. (13), (14), which it does, we have chosen not 
to do so. 

For completeness, Fig. 3 shows the development in the standard scenario at a later time (roughly 15 times the 
rest-frame lifetime - this is still an intermediate time due to the two gamma factors) - and Fig. 4 shows the standard 
scenario at an even later time (roughly 1400 times the rest frame life time), where equilibrium is almost reached. 

One should note that distribution functions are defined according to the cylinder coordinates used in the code - 
sec the Appendix VII [especially Eq. 24 which shows that the distribution function has dimension of time] for details. 
Also, one should note that p y in the plots are, in fact, the transverse momentum - not the momentum in one of the 
transverse directions. The unit of p in the figures is dp = ^ = which is the distance between adjacent bins in 
momentum space. 
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VI. CONCLUSION 



We have investigated a neutrino-pseudoscalar gas with only decays v\ — > is 2 + <p an d inverse decays to obtain 
equilibrium. We started with an anisotropic distribution of V\ and confirmed that kinetic equilibrium is delayed 
by two Lorentz 7-factors - one for time dilation of the heavy neutrino lifetime and one from the opening angle ie. 
from the transformation of the isotropic distribution of the decay products in the rest frame of the decaying particle 
back to the (cosmic) laboratory frame. We found explicit analytical expressions for the rate of creation of transverse 
momentum - both in a case with no background of the decay products and in the case of thermal backgrounds and 
the ultrarelativistic limits hereof. We have confirmed this behavior in numerical simulations as well - though we had 
to make a thermal smear of the initial anisotropy, making the analytical and numerical results open to a qualitative, 
but not quantitative, comparison. 

Furthermore we have obtained updated bounds on the neutrino-pseudoscalar coupling constant as well as on the 
neutrino lifetime in the realistic case of a thermal background of neutrinos and pseudoscalars. 



We would like to thank Steen Hannestad and Georg Raffclt for comments and discussion during the coding and 
writing part of this project. AB is supported by the Danish Council for Independent Research — Natural Sciences. 



In order to follow this numerically we notice that we have three particles in three momentum coordinates - that is, 
nine dimensions (no isotropy). We notice that when we assume Maxwcll-Boltzmann statistics particles 2 and ip behave 
alike. This means that if starting conditions are the same, we have to track only one of them. Even though there is 
not isotropy, azimuthal angles are arbitrary. So we end up with four dimensions, two particles, with a momentum in 
the initial beam direction and momentum in a transverse direction. So we want to integrate the remaining coordinates 
out. We start with 
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VII. APPENDIX: NUMERICS 



(2n f5 4 ( Pl -p2 -p v )g 2 m 2 {h - / 2 /„) 



2E 1 J (2n) 3 2E 2 (2tt) 3 2E v 

dp V xdp ip Rdp 2x dp2RP2RP v >R. 

(2tt) 2 E 1 E 2 E v , 



S(E 1 -p 2 -p v )5{pi x -P2x -P<px) 




(21) 



where we have defined (including aligning the coordinate system with particle 1) momenta 
Pi = {Pix,Pm,0),P2 = (p2x,P2RCOs(ip) 7 p 2 Rsm(ip)) ,p v = {p V x,P V RCOs(0) 7 p ipR sin(e)). 
After integration we find 



C[l} = 



dp vx dp v Rdp2xdp2RP2RP V R 

r K 2 E\E2E ip 



■5{Ei -P2- P v )5(pix - P2x - P V x)g 2 m 2 
Mr 1 



(h - hU) 



(22) 



where S is given by 



S = ±p\rP%r - (pIr - P%r - p\r) ■ 



For numerical purposes, let us underline the formula in the way it should be implemented. 



(23) 



however the vectors will not be introduced. Rather we use 



/ f(pl)dp0 = / f{Ptx,PtR,Pe)PiRdpe = 2irp lR * f(p lx ,p lR ,pe) = 27rp lR f(p lx ,p lR ) (24) 
J e Je 

and likewise for C[i\. This means that the function / that we implement is of dimension E^ 1 and its derivative 
dimensionless. The equation implemented is thus 

df(pix,PiR) /V M1 dp lpx dp lpR dp2xdp2RP2RP<pR xf „ s 
= 6 W = JF TP TP P V ) 



dt 1 1 w 3 E 1 E 2 E lfi 

8(pix -P2x-P v x)9 2 m 2 (f 1 - / 2 / v )^= 
1 

The implementation of 2 is quite easy since for fixed momenta pi,P2,P^> 



or 



(25) 



C[l] = -C[2] (26) 



27rp 1J jC[l] - -27rp2flC[2]. (27) 
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